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Abstract. We describe a Continuous Hugoniot Method for the efficient simulation of shock wave fronts. This approach 
achieves significantly improved efficiency when the generation of a tightly spaced collection of individual steady-state shock 
front states is desired, and allows for the study of shocks as a function of a continuous shock strength parameter, v p . This is, to 
our knowledge, the first attempt to map the Hugoniot continuously. We apply the method to shock waves in Lennard-Jonesium 
along the < 100> direction. We obtain very good agreement with prior simulations, as well as our own benchmark comparison 
runs. 



INTRODUCTION 

On-going experiments the University of Texas at Austin 
are investigating the shock-strength dependence of melt 
time scales in tin using terawatt laser systems on sub- 
picosecond time scales. Large numbers of trials explor- 
ing many different shock velocities will be possible. We 
present here a method intended to provide direct compar- 
isons with these experiments. 

There are two major difficulties with applying tradi- 
tional simulation methods [1] [2] [3] to the study of the 
dynamics near the shock front: (1) To produce the en- 
vironment at the front, one must simulate a large and 
ever-growing system, of which the front constitutes only 
a very small fraction; and (2) The conditions within a 
steady-state shock take long times to arise, and each 
computationally-expensive shock run results in only a 
single data point. 

The constrained dynamics methods of the Hugoniostat 
and others [4] [5] offer a solution to the first issue, but 
provide no information about non-equilibrium dynamics 
at the shock front that would be needed to compare with 
experiments. The approach of Zhakhovskii et al. [6] suc- 
ceeds in addressing the first point at the shock front, but 
does not address the second. We generalize and expand 
on their methods. First, we concentrate our efforts on the 
neighborhood of the shock interface, thereby increasing 
computational efficiency, and second, we map system re- 
sponse to a continuum of shock strength final states in a 
single run. These, combined, constitute the Continuous 
Hugoniot Method. 



CONTINUOUS HUGONIOT METHOD 

Our simulation method makes the Hugoniot the thermo- 
dynamic path of our simulation system. This has not been 
possible experimentally. Figure 1 contrasts the experi- 
mental loading paths (Rayleigh lines) to each state of the 
Hugoniot and the loading path which we will use to reach 
the same state points. 
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FIGURE 1. (left) The Hugoniot, as a collection of final 
shock states; and (right) as the state-to-state path of the Con- 
tinuous Hugoniot Method. 

We outline our simulation method in four stages: 
Benchmark - We begin each simulation with a tradi- 
tional shock wave computation. This initial long dura- 
tion, full-system run is allowed to continue to its final 
shock steady state. This steady state is our first point on 
the shock Hugoniot and serves to seed our subsequent 
computations. In Figure 1 this is represented by the lower 
Rayleigh line from the initial state to the beginning of the 
Hugoniot ramp. 



Reduce - We truncate the system to a fixed width neigh- 
borhood around the front by removing particles which 
are beyond a set distance behind the shock front. We 
impose a warm impactor (piston) boundary condition 
and match the thermodynamic statistics at the rear purge 
point using a strong Langevin thermostat. A buffer of 
pristine material is preserved ahead. The length of the re- 
duced system is determined by the system's thermaliza- 
tion length. The piston driving velocity and mean veloc- 
ity of the thermostat are equal, and this value, v p , serves 
as the external control parameter for shock strength. 
Shock strength ramp - To this reduced system, we intro- 
duce a quasistatic increase of the shock forcing parame- 
ter. The goal is to increase the shock strength, while al- 
ways maintaining a direct mechanical coupling between 
the forcing at the piston and the response at the shock in- 
terface. The piston velocity increases by a set amount per 
timestep. The temperature of the stochastic forcing is up- 
dated every purge to match the distribution at the purge 
point. The point of forcing always remains a set distance 
behind the shock front. This process continues until the 
shock strength parameter has reached its terminal value. 
Terminal benchmark comparison - Finally, a second 
benchmark run is made with a shock strength matched 
to the final value of the shock strength ramp. This allows 
the final state of the Hugoniot ramp to be compared di- 
rectly to a traditional shock run to identify any problems 
and to serve as an error check. 

Validity of Shock States 

The validity of our purge technique assumes that the 
back of our system is in equilibrium. We verify this by 
velocity distribution analysis before we reduce the sys- 
tem. We assume our shocks are always in steady state. 
This is guaranteed, so long as our ramp loading is qua- 
sistatic. If so, the small changes have time to equilibrate 
across the entire system and the driving and response are 
mechanically coupled. If not, then the foundation of the 
Hugoniot-Rankine equations is eroded, and off-Hugoniot 
states are produced. 

The essential property of this method is our abil- 
ity to very slowly ramp the shock strength parameter 
while maintaining correct thermodynamic conditions at 
a roughly constant distance behind the front. This is not 
currently achievable in experiment. Instead, in experi- 
ments driving forces are imposed at ever-increasing dis- 
tances. If one ramps the driving velocity in an experimen- 
tal situation, the result is isentropic compression rather 
than a shock response. 

Velocity Ramp Rate 

We can estimate an upper bound for the quasistatic ramp 
rate of the shock strength control parameter. The velocity 
is scaled by units of the wave speed of the compressed 
material, Cs, and the time is scaled by the return time of 



the acoustic waves in the system, 2L/Cs- We, thus, get 
the nondimensional condition for a quasistatic ramp. 



dv p _ d(v p /C s ) 



dt 



<1 



t<4 «> 



Note that the upper bound for quasistatic ramp rate goes 
to zero for large systems. The advantage of the method is 
lost when systems are too large. 



APPLICATION TO LENNARD-JONES 

As a prelude to more realistic but computationally inten- 
sive studies, we test these ideas with the Lennard-Jones 
potential, which has a well-documented solid shock re- 
sponse [7] [3]. 

Simulation Details 

We use the cubic-spline Lennard-Jones 6-12 potential 
[8] in order to allow easy comparison with published 
Hugoniot results of Germann et al. [1]. The shock was 
oriented along the <100> direction of the fee crystal 
with unit cell dimension 5.314 A = 1.561(7. Initial tem- 
perature was varied with a weak Langevin thermostat 
from zero to 10K = 0.083 k B T/e. Results are found not 
to depend on the initial temperature for shock driving ve- 
locities, v p above 0.75 C a . Systems were 20 x 20 lattice 
planes in cross-section with transverse periodic boundary 
conditions. The timestep was 0.3 femtoseconds. 

The traditional shock simulation runs (benchmark 
runs) were driven by a warm impactor and reached 
600 A in length (~ 100,000 particles). Continuous Hugo- 
niot Method runs were held to 200 A in length (~ 20,000 
particles), at any one time. The cumulative distance cov- 
ered by these treadmilling runs was almost 1.3 fim in 
length, and would have required over 1,200,000 particles 
in a conventional simulation. Every shock was given time 
to establish a steady state (usually 30 to 60 ps). The ramp 
rate of 0.001 m/s per step = 3.3 x 10 12 m/s 2 is used for 
the remainder of this article. We report on our efforts in 
the strong shock regime, for driving velocities ranging 
from v p = 0.75C o to 1.5C . 

Principal Hugoniot Results 

The Continuous Hugoniot Method allows a system to 
move directly from one shock state to another. Therefore, 
the path of our Continuous Hugoniot Method through 
U s -v p space during a single run is the principal Hugoniot 
of final shock states in the material. This is true to the 
extent that there is quick convergence within the reduced 
system to the values far behind the shock. Figure 2 shows 
such a path for a simulation which runs from v p = 0.75 C 
continuously through v p = 1.5C„. Initial and terminal 
benchmark runs, bookend the ramp. The Hugoniot fit 



proposed by Germann et al. is plotted in its applicable 
range. 



4 





1 1 1 1 1 1 1 1 1 1 1 


, | .... | , 

o 




— Simulation result 






— LJ HuQoniot fit (Germann, et al.) 






Full "benchmark" runs 
























... i .... i . . 


. 1 .... 1 . 



0.5 0.75 1 1.25 1.5 



particle velocity (Vp / C Q ) 
i 1 1 1 1 1 r 



1.6 - 




FIGURE 2. (top) Results of the Continuous Hugoniot 
Method, the fit of Germann et al. [1], and bookend bench- 
mark runs, (bottom) Density profiles snapshots from a shock 
strength ramp (solid), from the initial and terminal benchmarks 
(dashed), and from the Hugoniot fit (dotted). 

We see excellent agreement of our method's results 
with both comparisons. In the lower range of v p , where 
we can compare to the published fit, our data overlays 
the fit very nicely and continue it smoothly beyond the 
range for which it was originally published. At higher 
shock strengths our data stiffens (as it should), showing 
a super-linear increase in Us vs v p . In the upper range of 
v p , we compare to our terminal benchmark results, which 
also agree well. 

Comparison of Density Profiles 

Figure 2 shows a series of four density profile snap- 
shots taken from a continuous Hugoniot ramp (shown as 
solid lines). They are at v p — 0.75 C„, 0.9 C„, 1.2 C and 
1.5 C . The density profiles of the two benchmark runs 
and the final densities predicted by the Hugoniot fit for 
0.75 C and 0.9 C () are plotted. 

The profiles produced by the Continuous Hugoniot 
Method agree well with the results of the benchmark runs 
up to the point where they are purged. The benchmark 
densities, however, continue to grow beyond this point. 
In both cases the density predicted by the published 



fit is approximately 2% larger than the average final 
density produced by the Continuous Hugoniot Method. 
It appears that the density requires a larger spatial region 
than we have provided in order to converge completely 
to its asymptotic value. 

Comparison of Final States 

The final state of the Hugoniot ramp is a particularly im- 
portant point of comparison because it is the state of the 
maximum integrated error. Figure 3 provides snapshots 
from the terminal benchmark (top) and the Continuous 
Hugoniot Method (bottom). The slices are along the x-z 
plane at the final piston velocity v„ = 1.5C U . Both sys- 
tems exhibit a strong disordering transition on similar 
length scales, islands of incomplete disordering, and a 
common sharp density rise. Both have also developed 
forward-reaching features ahead of the front. Figure 3, 
shows very good agreement between the radial distribu- 
tion functions for the material behind each front. 
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FIGURE 3. Structural dynamics and rdf comparison - 

Particle slices in the region of the shock front for states from 
(middle) the Continuous Hugoniot Method and (top) the termi- 
nal benchmark. Also shown (bottom) is a comparison of radial 
distribution functions for the compressed material extending 10 
nm behind each front. v p = 1.5 C a . 

Efficiency and speed up 

The computational speedup provided by the Continuous 
Hugoniot Method depends upon the density of points 
with which one wants to locate the Hugoniot. We find 
that the computation time to compute two benchmark 
runs by traditional methods is approximately equal to 



the computation time we employed to compute the entire 
intermediate Hugoniot, via our method. Thus if we had 
chosen to trace out the Hugoniot by interpolation with 
10 conventional shock computations, we would have re- 
quired 5 times more computation. More generally, let 
Nh be the number of traditional runs needed to map 
the Hugoniot between two states as a function of shock 
strength. Then the speed up is linear in Nh, given roughly 
by Nh/2. Alternatively, we point out we were able to sim- 
ulate the cumulative effect of 1,280,000 particles with the 
resources required to hold only 20,000 at any one time. 



CONCLUSIONS 

We have presented the Continuous Hugoniot Method for 
efficient simulation of the dynamics in steady-state shock 
fronts. We confirmed in a Lennard-Jones system that the 
loading path of the Continuous Hugoniot Method fol- 
lowed the published Hugoniot fit. Comparison of parti- 
cle snapshots and radial distribution functions at the fi- 
nal state of the Continuous Hugoniot Method ramp also 
showed good agreement with traditional shock methods. 

All these measurements were made with greatly re- 
duced computational expenditure over traditional meth- 
ods. These savings are proving critical as the method is 
applied to more realistic and computationally costly po- 
tentials such as tin. 
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